Efficient Global Aerodynamic Modeling from Flight Data 

Eugene A. Morelli 1 

NASA Langley Research Center, Hampton, Virginia, 23681 


A method for identifying global aerodynamic models from flight data in an efficient 
manner is explained and demonstrated. A novel experiment design technique was used to 
obtain dynamic flight data over a range of flight conditions with a single flight maneuver. 
Multivariate polynomials and polynomial splines were used with orthogonalization 
techniques and statistical modeling metrics to synthesize global nonlinear aerodynamic 
models directly and completely from flight data alone. Simulation data and flight data from 
a subscale twin-engine jet transport aircraft were used to demonstrate the techniques. 
Results showed that global multivariate nonlinear aerodynamic dependencies could be 
accurately identified using flight data from a single maneuver. Flight-derived global 
aerodynamic model structures, model parameter estimates, and associated uncertainties 
were provided for all six nondimensional force and moment coefficients for the test aircraft. 
These models were combined with a propulsion model identified from engine ground test 
data to produce a high-fidelity nonlinear flight simulation very efficiently. Prediction testing 
using a multi-axis maneuver showed that the identified global model accurately predicted 
aircraft responses. 


Nomenclature 


a x ,a y ,a z 

- 

body-axis translational accelerometer measurements, g 

AirSTAR 

= 

Airborne Subscale Transport Aircraft Research 

b 

= 

wing span, ft 

c 

= 

wing mean aerodynamic chord, ft 

c x ,c Y ,c z 

= 

body-axis nondimensional aerodynamic force coefficients 

Q' Cf 1 

= 

body-axis nondimensional aerodynamic moment coefficients 

GPS 

= 

global positioning system 

INS 

= 

inertial navigation system 

> I y > ’ I xz 

= 

mass moments of inertia, slug-ft 2 

m 

= 

aircraft mass, slug 

M t 

= 

body-axis pitching moment from engine thrust, ft-lbf 

p, q,r 

= 

body-axis roll, pitch, and yaw rates, rad/s or deg/s 


= 

dynamic pressure, lbf/ft 2 

S 

= 

wing reference area, ft 2 

SIDPAC 

= 

System IDentification Programs for AirCraft 

T 

= 

maneuver length, sec 

T T 

*x> ± z 

= 

body-axis engine thrust, lbf 

V 

= 

true airspeed, ft/s 

a 

= 

angle of attack, rad or deg 

P 

= 

sideslip angle, rad or deg 

S e ,S a ,S r 

= 

elevator, aileron, and rudder deflections, rad or deg 

<j), 0,y/ 

= 

Euler roll, pitch, and yaw angles, rad or deg 

I 

= 

covariance matrix 
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superscripts 

T 

-1 

subscripts 

eg 

o 


transpose 
estimate 
time derivative 
matrix inverse 


center of gravity 
reference value or base term 


I. Introduction 

A CCURATE flight simulation has many important applications in aircraft dynamics and control, such as pilot 
training, mission rehearsal, dynamic analysis, and control system design 1 ' 3 . Typically, flight simulations are 
based on wind tunnel data and/or results from computational methods 2 ' 8 . Flight data are often used subsequently to 
modify or augment simulations for improved fidelity to the real flight vehicle 9 ' 12 . Furthermore, the Federal Aviation 
Administration (FAA) requires flight data matching criteria to be satisfied as part of the certification process for 
commercial flight training simulators. 

Resources are always limited, and wind tunnel testing and aerodynamic calculations require time and money. 
Building an aircraft simulation based only on flight data, without wind tunnel testing or aerodynamic calculations, 
could be effective for certain types of rapid, relatively low-budget research and development projects, such as 
subscale flying models and unmanned aerial vehicles (UAV). In such cases, there is no human life risked aboard the 
aircraft, and typically the aircraft are known to be fly able. The capability to build a high-fidelity nonlinear 
simulation based on a few test flights would be very beneficial in terms of project cost and schedule. 

Using flight data to update a simulation database generated from wind tunnel data or aerodynamic calculations is 
often a time-consuming and iterative task done in an ad hoc way, although recent work has focused on automated 
simulation updating methods 11,12 . If an accurate global aerodynamic model could be generated from flight data 
alone, the problem of flight-updating a simulation database based on wind tunnel data or aerodynamic calculations 
could be avoided altogether. Even in cases where a simulation has been generated based on wind tunnel data and 
aerodynamic calculations, a flight-determined global aerodynamic model could be used for simulation validation, 
updating based on flight data, and quantifying the effects of Reynolds number and geometry differences between the 
full-scale flight vehicle and the wind tunnel test article. 

In this work, a method for identifying a global aerodynamic model efficiently from flight data alone is examined. 
The method includes a novel approach to designing the maneuvers used to collect flight data for modeling purposes, 
and a novel approach to global aerodynamic modeling based on the flight data. The next section describes the 
maneuver design, which leverages past work on input design for real-time dynamic modeling 13 ' 17 . Following this, 
the global aerodynamic modeling approach is described. The approach combines multivariate polynomials, 
polynomial splines, orthogonal function modeling theory, and statistical modeling metrics to identify global 
aerodynamic models for nondimensional force and moment coefficients. The novel flight test maneuver design and 
global aerodynamic modeling method are then demonstrated using simulated data and flight data from a subscale 
twin-engine jet transport aircraft. A flight simulation was created by combining previously developed simulation 
software 6,14 written in MATFAB® with the global aerodynamic model identified from flight data and a propulsion 
model identified from ground test data. A multi-axis doublet sequence maneuver executed at a single flight 
condition (nominal angle of attack) was used for prediction testing to evaluate the fidelity of the global aerodynamic 
models identified from flight data. Details of the flight test and prediction results are presented in Section VI. 

The T-2 subscale jet transport aircraft used for this study, described in detail in Section IV, was also tested 
extensively in the wind tunnel 5,8 , although the wind tunnel test article differed in some geometric details from the 
aircraft used for the flight tests. The wind tunnel data provided an aerodynamic database for the nonlinear 
simulation, which was used to validate the effectiveness of the global aerodynamic modeling procedure using flight 
data. Comparisons using global aerodynamic models identified from flight data and the wind tunnel aerodynamic 
database are shown in Sections V and VI. 

All of the experiment design, data analysis, and modeling tasks included in this work were done using system 
identification software written in MATFAB®, called System IDentification Programs for AirCraft, or SIDPAC 14 . 
SIDPAC is bundled with Ref. [14], and is therefore publicly available. The SIDPAC software toolbox was 
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developed at NASA Langley, and is continually expanded and improved. SIDPAC has been applied successfully to 
a wide variety of flight and wind tunnel experiments at NASA Langley 14 and elsewhere, and is used at more than 80 
institutions worldwide 18 . 


II. Maneuver Design 

A. Orthogonal Optimized Multi-Sine Input Design 

This section describes how orthogonal optimized multi-sine inputs are designed and why this particular input 
form is efficient for identifying dynamic models from flight data. More details on this input design technique and 
flight applications can be found in Refs. [13] -[17]. 

The general idea is to excite the aircraft using perturbation inputs with wideband frequency content over a range 
of frequencies that encompasses the expected modal frequencies for the aircraft dynamic response. The excitations 
are implemented as perturbations to the control surface deflections by summing designed perturbation inputs with 
the actuator commands from the pilot and feedback control system, just before the actuator limiting on command 
rate and position. 

Each designed perturbation input is a sum of sinusoids with unique frequencies, optimized phase shifts, and 
specified power distribution. Component frequencies are selected to cover a frequency band of interest, similar to 
frequency sweeps. The wide -band frequency content of the inputs is important because there is naturally some 
uncertainty as to what the modal frequencies are for the aircraft in flight. Wide -band inputs provide robustness to 
this uncertainty. Phase shifts for the sinusoidal components of each input are optimized to achieve low peak-to-peak 
amplitude and high input energy content for the sum of sinusoids. Amplitudes of the individual sinusoidal 
components can be chosen to achieve a specific power distribution. 

Multiple inputs are designed to be mutually orthogonal in both the time domain and the frequency domain, and 
are optimized for maximum data information content in multiple axes over a short time period, while minimizing 
excursions from the nominal flight condition. The mutual orthogonality of the inputs allows simultaneous 
application of multiple inputs, which helps to minimize excitation time, but more importantly for this work, provides 
continuous multi-axis excitation as the aircraft flies through time-varying or precarious flight conditions. 

Each perturbation input Uj , which is to be applied to the / h individual control surface, is comprised of a set of 

summed harmonic sinusoids with individual phase shifts (j) k , 


u i= Z A k sin \—f- + 4k} j = l,2,...,m (1) 

where M is the total number of available harmonically -related frequencies, T is the time length of the excitation, and 
A k is the amplitude for the k th sinusoidal component. The variable t represents a vector of N discrete time 

points, t - [f (0) t(l ) ... /(A -l)] 7 , and Uj represents the vector of corresponding amplitudes for the / h input, 

u j = | ~Uj (0) Uj (l) ... Uj (A-l)] . Each of the m inputs is comprised of selected components from the pool 
of M harmonic sinusoids with frequencies o\ -l^k/T , k = 1,2,...,M , where co M represents the upper 

limit of the frequency band for the excitation. The interval \ct)\ f co M ] rad/s specifies the range of frequencies where 
the aircraft dynamics are expected to lie. 

If the phase angles (/) k in Eq. (1) were chosen at random on the interval (-tc,tz\ rad, then in general, the various 
harmonic components would add together at some points to produce an input Uj with relatively large amplitude 

excursions. This is undesirable, because it can result in the dynamic system being moved too far from the reference 
condition selected for the experiment. To prevent this, the phase angles (j) k for the selected harmonic components 
are chosen to minimize relative peak factor RPF 14 , defined by 
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Relative peak factor is a measure of the efficiency of an input for dynamic modeling purposes, in terms of the 
amplitude range of the input divided by a measure of the input energy. The relative peak factor is scaled so that any 
individual sinusoidal component (such as any one of the summands in Eq. (1)) has RPF = 1. Low relative peak 
factors are desirable and efficient for estimating dynamic model parameters, because the objective is to excite the 
dynamic system with good input energy over a variety of frequencies while minimizing the input amplitudes in the 
time domain, to avoid driving the dynamic system too far away from the reference condition. 

For a composite signal with more than one sinusoidal component, as in Eq. (1), the goal of designing an input 
with minimum RPF is achieved by adjusting the phase parameters </> k for the sinusoidal components of the input. 
The resulting optimization problem is non-convex; however, a simplex algorithm 14,19 can be applied to find a 
solution. 

The integers k specifying the frequencies for the / h input Uj are selected to be unique to that input, but are not 

necessarily consecutive. A good approach for multiple inputs is to assign integers k to each input alternately. This 
is illustrated in Figure 1 for a flight test maneuver design on the T-2 subscale jet transport aircraft described in 
Section IV. In that case, there were 3 inputs: elevator, rudder, and aileron, and a total of 30 frequencies (M = 30) . 

The frequencies were interleaved among the three inputs to achieve wide -band frequency content for each input. 
This provided robustness to uncertainty in how each control excites the dynamic modes of the aircraft. Because 
each input has wide-band frequency content, the same input design can be applied at various flight conditions, which 
simplifies the flight test and reduces flight computer memory requirements. It is even possible to use the same input 
design for different aircraft, because of the wide-band frequency content of the excitation inputs. 

To achieve a uniform power distribution, the A k are selected as 

\=A= Vk (3) 

>Jn 

where n is the number of sinusoidal components included in the summation of Eq. (1) for i# • , and A is the 
amplitude of the composite input Uj . Therefore, with uniform power distribution, selection of the A k reduces to 

selecting a single value for the input amplitude A . Each 
input Uj can of course have arbitrary amplitude A, 

subject to practical flight testing and modeling 
constraints. 

It is also possible to modify the power at individual 
frequencies for each input, to focus the excitation on 
frequencies near where the natural frequencies of the 
dynamic modes are believed to be, or to avoid exciting 
structural responses, for example. For each input, the 
power spectrum can be tailored by selecting the A k in 
Eq. (1) to distribute power over the spectral components. 

The power spectra shown in Figure 1 are normalized, so 
the effects of individual control surface amplitudes are 
excluded. This means that for each input, the sum of all 
the spectral line ordinates (sum of the heights of the bars 
for each input) is 1 . 

When the frequency indices k selected for each 
input Uj in Eq. (1) are distinct from those chosen for the 

other inputs, then the frequency content of each Uj consists of distinct spectral lines in the frequency domain, as can 
be seen in Figure 1 . Therefore, the vectors of Fourier transforms for the inputs as a function of frequency have inner 



frequency (Hz) 

Figure 1 . Multiple orthogonal phase-optimized 
multi-sine input spectra 
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products equal to zero. In this sense, the inputs are mutually orthogonal in the frequency domain, because each 
input contains frequencies that are not present in the other inputs. 

In the time domain, a sum of harmonic sinusoids is orthogonal to any other sum of sinusoids with harmonically- 
related frequencies, regardless of the constant phase shift of each sinusoidal component 14,15 . Consequently, the 
inputs are also mutually orthogonal in the time domain. 

An objective for the experiment design is to excite the aircraft dynamics in all axes over a short time period by 
moving multiple control surfaces simultaneously. This is particularly important in situations where the reference 
flight condition is time- varying, or cannot be maintained for very long. 

Since more than one surface is being moved, it is advantageous for modeling purposes if the Uj applied to the 

control surfaces are mutually orthogonal. This helps the dynamic modeling by completely de -correlating the inputs, 

which improves the accuracy of control effectiveness 
estimates. Using the input design method described here, it 
is possible to make all of the Uj mutually orthogonal in 

both the time and frequency domains, while also 
minimizing relative peak factor for each Uj , which keeps 

the aircraft from departing significantly from the reference 
flight condition. This gives the analyst the flexibility to use 
either time-domain or frequency-domain parameter 
estimation methods, while retaining the desirable feature of 
mutually orthogonal inputs. 

Figure 2 shows the perturbation input time series for a 
maneuver design flown on the T-2 subscale jet transport 
aircraft. These inputs are mutually orthogonal in both the 
time and frequency domains. The inputs were computed 
from Eq. (1) and the information in Table 1, where the 
phase angles fa were optimized for minimum relative 
peak factor. Because of the various frequencies and phase 
angles, and the small amplitudes of the perturbation inputs, 
applying these inputs simultaneously to the aircraft 
produces a dynamic response similar to what might be seen 
in flight through light to moderate turbulence. The aircraft 
stays near the reference condition, but responds 
dynamically about that condition. In practice, pilot inputs and feedback control can act to ruin the input 
orthogonality; however, good modeling results require only low correlations, not zero correlations, so that slightly 
imperfect inputs still work quite well. 

B. Maneuver Design for Efficient Global Aerodynamic Modeling 

In past work 13 ' 17 , inputs such as those shown in Figure 2 were applied at a selected flight condition. The pilot 
flew the aircraft to the reference flight condition, then activated an automated excitation system that added inputs 
like those shown in Figure 2 to the trim control deflections coming from the pilot and feedback control system. This 
produced excellent data for aerodynamic parameter estimation at the selected flight condition. This approach could 
be called local aerodynamic modeling. 

For global aerodynamic modeling, one conventional approach is to combine local aerodynamic modeling results 
obtained from local perturbation maneuvers. This requires numerous and accurate acquisitions of particular flight 
conditions, followed by the application of perturbation excitations at each condition, and combining the local results 
to produce a global model. 

More efficient global aerodynamic modeling can be achieved by continuously applying multi-axis perturbation 
inputs while the aircraft flight condition is varied slowly. This approach is practical because the orthogonal 
optimized multi-axis perturbations excite the aircraft dynamics in a very time -efficient manner with high data 
information content, so that the aircraft dynamics can be sufficiently excited even when the flight condition is 
changing. 

To implement an efficient global aerodynamic modeling maneuver, the pilot began with a steady wings -level 
trim condition at low angle of attack, then initiated automated multi-axis excitation inputs like those shown in 




time (s) 

Figure 2. Multiple orthogonal phase-optimized 
multi-sine excitation inputs 
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Figure 2. While the excitation inputs were being applied, the pilot pulled back slowly on the stick to induce a slow 
increase in the nominal angle of attack. The multi-axis excitation inputs were continuously applied additively to the 
control surface deflections commanded by the pilot and control system. Pilot inputs on lateral stick and rudder pedal 
were essentially zero, but the elevator deflection commanded by the pilot changed slowly to implement the slow 
increase in nominal angle of attack. An example of the resulting flight data is shown in Figure 3. This maneuver 
produces very informative data over a wide range of nominal angle of attack. Only subsonic aerodynamics at 
relatively low altitude were being studied, so the effects of changing airspeed were adequately modeled by 
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Figure 3. T-2 flight data for a global aerodynamic modeling maneuver 


conventional nondimensionalization using dynamic pressure. Since the nominal angle of attack changed slowly, this 
maneuver can be considered a combination of informative multi-axis excitation data for many different nominal 
angles of attack, executed in a single, efficient, combined maneuver. 

Figure 4 shows cross-plots of aircraft states and controls using data from the maneuver shown in Figure 3. 
These plots demonstrate that a wide range of the explanatory variables generally used for aerodynamic modeling 
was swept through during this single maneuver. Note also that the cross -plots generally do not show diagonal lines 
or ellipses, which means that the explanatory variable data from this maneuver had very low pair-wise correlations. 
Low pair-wise correlations mean that the aerodynamic dependencies on the explanatory variables can be identified 
accurately and without ambiguity. Cross-plots for other aircraft states and controls were similar in that the plots 
indicated low pair-wise correlations for the explanatory variables. 
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Figure 4. T-2 flight data cross plots for a global aerodynamic modeling maneuver 


This novel maneuver exhibits a combination of low pair-wise correlations, multi-axis excitation, and slowly- 
varying flight conditions that cover a large portion of the explanatory variable subspace for aerodynamic modeling. 
These characteristics make the maneuver very effective and efficient for global aerodynamic modeling. 
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A series of these maneuvers were flown on the T-2 subscale jet transport aircraft. Each maneuver had a constant 
power setting and aircraft configuration, and all used multi-axis perturbation time series, similar to what is shown in 
Figure 2. The wide-band frequency content of the excitation inputs made them effective for a wide variety of flight 
conditions, power settings, and aircraft configurations. Over several flight test deployments, some variations in the 
frequency content, time length, and amplitudes of the excitation inputs were implemented, in order to study the 
effects of these variations from an experiment design standpoint. Table 2 lists the power setting and aircraft 
configuration for each maneuver in the series flown on the T-2 aircraft. 

III. Global Aerodynamic Modeling 

The objective of global aerodynamic modeling is to identify a model for each nondimensional aerodynamic force 
and moment coefficient as a function of explanatory variables that can be measured, such as angle of attack, pitch 
rate, and control surface deflections, over a large range of the explanatory variables. There are two main difficulties 
in doing this: 1) experiment design to collect dynamic data over a large range of the explanatory variables, which 
was addressed in Section II, and 2) identifying an accurate global model, which is the subject of this section. 

In the typical approach to global aerodynamic modeling, a series of linear or simplified local models, which are 
valid for relatively small ranges of the explanatory variables, are joined together in some way to implement a global 
aerodynamic model. In another approach, a multivariate orthogonal function modeling technique was developed to 
identify global models, which were generally nonlinear and valid over relatively large ranges of the explanatory 
variables 20 . However, for very large ranges of the explanatory variables, or severe local nonlinearities, this approach 
in some cases compromises local model fit in order to achieve a better global model fit. This is the result of using 
global orthogonal polynomial functions for the modeling, which do not have the ability to change locally without 
modifying the entire model. Other approaches have used data partitioning with simplified local model structures 21 
and localized modeling functions 22 to address this problem, and also to make local model updates easier. However, 
these approaches require a procedure or analyst judgment for: 1) partitioning the explanatory variable space; 
2) selecting the location of the local modeling functions in the explanatory variable space; and/or 3) selecting the 
mathematical structure of the local models. 

In this work, the multivariate orthogonal function modeling first described in Ref. [20] was applied using both 
the explanatory variables and spline functions of the explanatory variables. This approach provides the required 
local nonlinear modeling capability, while retaining easy physical interpretation of the model, and automating the 
selection of optimal modeling complexity necessary to accurately characterize the functional dependencies. 

The next sections describe the approach to global aerodynamic modeling using multivariate orthogonal functions 
derived from multivariate polynomials and spline functions of the measured explanatory variable data. 

A. Aircraft Aerodynamic Modeling 

For global aerodynamic modeling from flight data, nondimensional aerodynamic force and moment coefficients 
are used as the response variable (also called the dependent variable) in the modeling problem. A separate modeling 
problem is solved for each force or moment coefficient, corresponding to minimizing the equation error in each 
individual equation of motion for the six rigid-body degrees of freedom of the aircraft. Values for the 
nondimensional aerodynamic force and moment coefficients cannot be measured directly in flight, but instead must 
be computed from measured and known quantities using the following equations 14 
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These expressions retain the full nonlinear dynamics in the aircraft equations of motion. 

The result is TV values of the nondimensional force and moment coefficients, where TV is the number of data 
points. These values are often called measured force and moment coefficients, even though they not measured 
directly, but rather computed from other measurements and known quantities. Explanatory variables such as angle 
of attack, Mach number, pitch rate, and control surface deflections, are measured directly. 

The desired form of the global aerodynamic model is a mathematical model structure with estimated model 
parameter values and associated uncertainties, relating the nondimensional aerodynamic force and moment 
coefficients to aircraft states and controls that can be measured. 


B. Multivariate Orthogonal Function Modeling 

The form of a multivariate orthogonal function model is 


z = a 1 p l +a 2 p 2 + ... + a n p n +e (9) 

where z is an TV-dimensional vector of the response variable (e.g., nondimensional force or moment coefficient), 
Z = [zi>Z2>—> z n] T * modeled in terms of a linear combination of n mutually orthogonal modeling functions 
p j, j = l,2,...,n . Each pj is an TV-dimensional vector which in general depends on the explanatory variables. The 
cij , j = 1,2,..., n are constant model parameters to be determined, and £ denotes the modeling error vector. 

Equation (9) represents a mathematical model used to represent functional dependencies in the measured data. 
The important questions of determining how the modeling functions pj are computed from the explanatory 

variables, as well as which modeling functions should be included in Eq. (9), which implicitly determines n , will be 
addressed later. At this point, the properties of a multivariate orthogonal function model are examined. 

Define an Nxn matrix P , 


P = [Pl>P2’-’Pn] ( 10 > 

and let a = \a l ,a 2 ,—,a n ^ . Equation (9) can then be written as a standard least squares regression problem, 

z = Pa + s (11) 


The error vector s is to be minimized in a least squares sense. The goal is to determine a that minimizes the least 
squares cost function 


J = —(z~ Pa) T (z-Pa)=—£ T £ 


( 12 ) 


The parameter vector estimate a that minimizes this cost function is computed as 1 


d=[p T pJ x 


P T Z 


(13) 


The estimated parameter covariance matrix is 1 




(d-a)(d-a) T =cr 2 ^P T P^ 1 


( 14 ) 
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( 15 ) 


where E is the expectation operator, and the fit error variance cj 1 can be estimated from the residuals 

v =z- Pa 


using 


1 

(N-n) 


(z-Paf \z-Pa) 


(N-n) 


(16) 


Parameter standard errors are computed as the square root of the diagonal elements of the 27^ matrix from Eq. (14), 

„ 9 

using cj from Eq. (16). The identified model output y is computed as 


y = Pa 


(17) 


In conventional least squares modeling, the modeling functions (columns of P) are often polynomials in the 
explanatory variables. This approach corresponds to using the terms of a multivariate Taylor series expansion to 
approximate the functional dependence of the response variable on the explanatory variables. If the modeling 
functions are instead multivariate orthogonal functions generated from the explanatory variable data, it is easier to 
determine an appropriate model structure, because the explanatory capability of each modeling function is 
completely distinct from all the others. This decouples the least squares modeling problem, as will be shown now. 

For mutually orthogonal modeling functions, 


pj Pj = 0 , i=tj , i, j = 1 , 2 ,..., n 


(18) 


and P T P is a diagonal matrix with the inner product of the orthogonal functions on the main diagonal. Using 
Eqs. (10) and (18) in Eq. (13), the / h element of the estimated parameter vector a is given by 


“j=(p r j z )l(p T i Pj) 


Using Eqs. (10), (18), and (19) in Eq. (12), 


J = 


1 

2 




(19) 


( 20 ) 


Equation (20) shows that when the modeling functions are orthogonal, the reduction in the least square cost 
function resulting from including the term aj pj in the model depends only on the response variable data z and the 

added orthogonal modeling function pj . The least squares modeling problem is therefore decoupled, which means 

each orthogonal modeling function can be evaluated independently in terms of its ability to reduce the least squares 
model fit to the data, regardless of which other orthogonal modeling functions are already selected for the model. 
When the modeling functions pj are instead polynomials in the explanatory variables (or any other non-orthogonal 

function set), the least squares problem is coupled, and iterative analysis is required to find the subset of modeling 
functions for an adequate model structure. 

The orthogonal modeling functions to be included in the model are chosen to minimize predicted squared error, 
PSE, defined by 23 


PSE = 


( z-Paf ( z-Pa ) 


N 


+ 


N 


( 21 ) 
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or 


PSE = — + a 2 max — 


( 22 ) 


The constant cr max is the upper-bound estimate of the squared error between future data and the model, i.e., the 
upper-bound mean squared error for prediction cases. The upper bound is used in the model over-fit penalty term to 
account for the fact that PSE is calculated when the model structure is not correct, i.e., during the model structure 
determination stage. Using the upper bound is conservative in the sense that model complexity will be minimized as 
a result of using an upper bound for this constant in the penalty term. Because of this, the value of PSE computed 
from Eq. (22) for a particular model structure tends to overestimate actual prediction errors on new data. Therefore, 
the PSE metric conservatively estimates the squared error for prediction cases. 

A simple estimate of cr^ax that is independent of the model structure can be obtained by considering <r^ x to be 
the residual variance estimate for a constant model equal to the mean of the measured response values, 


(i n 


i N 

= — Y\z i -z } 2 


(23) 


where 


z = 



(24) 


The PSE in Eq. (22) depends on the mean squared fit error, 2 j/N , and a term proportional to the number of 
terms in the model, n . The latter term prevents over-fitting the data with too many model terms, which is 
detrimental to model prediction accuracy 4,14,23 . While the mean squared fit error 2 j/N must decrease with the 


addition of each orthogonal modeling function to the model (by Eq. (20)), the over-fit penalty term <J^ nax n/N must 
increase with each added model term ( n increases). Introducing the orthogonal modeling functions into the model in 

order of most effective to least effective in reducing the mean squared fit error (quantified by [^p T jZ ) j{^p T j p for 

the / h orthogonal modeling function) means that the PSE metric will always have a single global minimum. 

Figure 5 depicts this graphically, using actual modeling results from Ref. [4] . The figure shows that after the 


first 6 modeling functions, the added model complexity 
associated with an additional orthogonal modeling 
function is not justified by the associated reduction in 
mean squared fit error. This point is marked by a 
minimum PSE , which defines an adequate model 
structure with good predictive capability. Ref. [23] 
contains further statistical arguments and analysis for 
the form of PSE given in Eq. (22), including 
justification for its use in modeling problems. 

Using orthogonal functions to model the response 
variable makes it possible to evaluate the merit of 
including each modeling function individually , using 
the predicted squared error PSE. The goal is to select a 
model structure with minimum PSE , and the PSE 
always has a single global minimum for orthogonal 
modeling functions. This makes the model structure 
determination a well-defined and straightforward 
process that can be (and was) automated. 


x 10 5 



Orthogonal Function Number 


Figure 5. Model structure determination using 
orthogonal functions and PSE 
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C. Generating Orthogonal Modeling Functions 

Multivariate orthogonal functions can be generated from ordinary multivariate functions in the explanatory 
variables using a Gram-Schmidt orthogonalization procedure. This approach is described in Refs. [4] and [14], 
which are the basis for the material presented here. 

The process begins by choosing one of the ordinary multivariate functions as the first orthogonal function. 
Typically, a vector of ones (associated with the bias term in the model) is chosen as the first orthogonal function, 

Pi= 1 (25) 

In general, any function of the explanatory variables can be chosen as the first orthogonal function, without any 
change in the procedure. To generate the next orthogonal function, an ordinary multivariate polynomial function is 
made orthogonal to the preceding orthogonal function(s). Define the / h orthogonal function pj as 


j - 1 

Pj = </ ~ X Yk} Pk J = 2 > 3 > •••> n t ( 26 ) 

k = 1 

where is the / h ordinary multivariate function vector. For example, each g j could be some ordinary polynomial 
function of the explanatory variables or a spline function of the explanatory variables. The y k j for k = 1, 2, ...,j — 1 

are scalars determined by multiplying both sides of Eq. (26) by pf , then invoking the mutual orthogonality of the 
p k , k = 1,2, ..., j , and solving for y k j 


Pk%i 

r k j= J y L k = l, 2,..., j-1 (27) 

PkPk 

The same process can be implemented in sequence for each ordinary multivariate function £ •, j - 2, 3, ..., n t . 

The total number of ordinary multivariate functions used as raw material for generating the multivariate orthogonal 
functions, including the bias term, is n t . It can be seen from Eqs. (25)-(27) that each orthogonal function can be 
expressed exactly in terms of a linear expansion of the original multivariate functions. The orthogonal functions are 
generated sequentially by orthogonalizing the original multivariate functions with respect to the orthogonal 
functions already computed, so that each orthogonal function can be considered an orthogonalized version of an 
original multivariate function. 

The multivariate orthogonal function generation method described here normally starts by generating all possible 
ordinary multivariate polynomials in the explanatory variables, up to a selected maximum order. For example, when 
modeling the nondimensional vertical force coefficient C z , the explanatory variables might be angle of attack a , 
nondimensional pitch rate q = qc/ 2V , and elevator deflection 5 e . If the selected maximum order is 3, then the 
ordinary multivariate polynomial modeling functions used as the raw material for the orthogonalization process 
would include terms like a,q,S e ,a ,S e ,aq ,aqS e , etc. Note that considering any other candidate explanatory 
variables, such as sideslip angle /? , can be done by simply including the sideslip angle among the group of 
explanatory variables. If it turns out that the sideslip angle is not needed to model C z , the model structure 
determination using orthogonal modeling functions will not select any orthogonal functions associated with sideslip 
angle. This occurs naturally and automatically in the course of the model structure determination process described 
earlier. Therefore, there is no harm in including explanatory variables that might not be important, except that 
additional computation time will be required to identify the model structure, because additional multivariate 
orthogonal functions will be generated and sorted. Similarly, if the maximum order is chosen higher than necessary, 
the only penalty would be the increased computation time necessary for generating and sorting the additional 
orthogonal functions. The final identified model would be the same. Consequently, the choices that the analyst 
needs to make are very easy and not critical to the quality of the final results. 
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If the p j vectors and the fy vectors are arranged as columns of matrices P and X , respectively, and the y k j 
are elements in the k th row and / h column of an upper triangular matrix G with ones on the diagonal, 



"1 

Yn 

ri3 * 

•• rin, 


0 

1 

Y 23 ' 

■■ Tin, 

G = 

0 

0 

1 - 

■■ Tin, 


(28) 


Then 


0 0 0 


X = PG 


which leads to 


P = XG~ l 


(29) 


(30) 


The columns of G 1 contain the coefficients for expansion of each column of P (i.e., each multivariate orthogonal 
function) in terms of an exact linear expansion in the original multivariate functions in the columns of X . 
Equation (30) can be used to express each multivariate orthogonal function in terms of the original multivariate 
functions. The manner in which the orthogonal functions are generated allows them to be decomposed without 
ambiguity into an expansion of the original multivariate functions, which have physical meaning. 

D. Conversion to Physically -Meaningful Multivariate Function Models 

After the model structure is determined using multivariate orthogonal modeling functions for minimum PSE , the 
identified model output is 


y=Pd (31) 

where the P matrix now includes only the n orthogonal functions selected in the model structure determination, 
n<n t . Each retained orthogonal modeling function can be decomposed without error into an expansion of the 

original multivariate functions in the explanatory variables, using the columns of G 1 in Eq. (28) corresponding to 
the retained orthogonal functions. Common terms are combined using double precision arithmetic to arrive finally 
at a model using only original multivariate functions in the explanatory variables. Terms that contribute less than 
0.1 percent of the final model root-mean- square magnitude are dropped. 

The final form of the model is a sum of ordinary multivariate functions in the explanatory variables, with 
associated model parameter estimates. Examples of the final model forms obtained are given later in the Results 
sections. 

E. Including Spline Functions in the Orthogonalization 

For global aerodynamic modeling, functional dependencies of the aerodynamic coefficients on the explanatory 
variables can exhibit significant localized variation. In these cases, a global polynomial model can be inadequate for 
capturing those local variations. To solve this problem, spline functions in the explanatory variables can be 
introduced as additional pseudo -explanatory variables. Splines have local modeling capability that global 
polynomial modeling functions do not have, because splines are polynomials defined only on selected intervals. 
Low-order polynomial terms defined on limited intervals can approximate nonlinearities quite well. Using 
polynomial splines also retains clear physical interpretation in the final model. 

Spline functions are defined as piecewise polynomials functions of degree m in one or more explanatory 
variables. The term “piecewise” means that the polynomial is different for specific ranges of the explanatory 
variables. Spline function values and derivatives agree at the points where the piecewise polynomials join. These 
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points are called knots, and are defined as specific values of each explanatory variable. A polynomial spline S m (x) 
of degree m with continuous derivatives up to degree m-1 , for a single explanatory variable v e [^O’-'Wx] > can 
expressed as 


where 


m k 

S m (*) = Z ' C r X T + D t (x ~ X; )“ 


r=\ i = 1 


(x-Xj) 



X > x t 
X<Xi 


(32) 


(33) 


and the C r and are constants and (x-x t )™ are the piecewise parts of the spline function. The values 
x lf x 2 , . . x k are knots which satisfy the condition 

X 0 <X l <X 2 <...<X k <Xma X (34) 


Three piecewise spline types are sketched in Figure 6 for the same selected knots in angle of attack. 

Note that when the spline knots are the same, a higher-order piecewise spline in a single explanatory variable can 
be computed as a multiplication of lower-order piecewise splines in the same explanatory variable, 

(* - Xi )™ = (x - Xj (x-Xif + m> 1 (35) 

This fact can be used advantageously in generating multivariate orthogonal modeling functions with excellent 
local modeling capability. This is best demonstrated with a simple example. 

Suppose that the nondimensional vertical aerodynamic force coefficient C z is being modeled with explanatory 
variables angle of attack a and elevator control deflection 8 e . For a maximum selected model order 2, the ordinary 
polynomial modeling functions that will serve as raw material for the orthogonalization process are: 

l, a , a 2 , aS e , 8 e , S 2 (36) 


which would lead to a final multivariate polynomial model of the form 


c z ~ c z n +c z r/ a + c z 9 OC 1 + CV a8 e +c z S e +C z S 2 


(37) 


where the values of the model parameters such as C z ° and C Zag would be estimated from the data using the 

orthogonal function modeling and subsequent decomposition procedures described earlier, and some of the terms 
might not be present, depending on the results of the model structure determination using orthogonal functions. This 
represents a global polynomial model. 

Now introduce a first-order spline term in angle of attack with a single knot located at 10 deg. Using the same 
maximum model order 2, the set of ordinary polynomial modeling functions expands to 

l,a,a 2 ,aS e ,S e ,Sj,(a-10) 1 + ,a(a-l0) 1 + ,S e (a-l0) 1 + ,(a- 10); (38) 


which would lead to a final multivariate polynomial model of the form 

c z = c z n + c z„ a + c z , q2 + Cz„, aS e + C z 8 e +c z S 2 


+ C Zi (a-10f + +C z x a(a-\Qi)\ + C z ] 8 e (a-lo£ + C Z 2 (a-10)J 


( 39 ) 
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As before, some of the terms in Eq. (39) might be not be present, depending on the results of the model structure 
determination using orthogonal functions. 

With this expanded set of ordinary polynomial 
modeling functions, more orthogonal functions will be 
generated, and there is now capability to accommodate 
changes in the linear dependence of C z on a for 

a > 10 deg , as well as additional a and aS e 
nonlinearities that only take effect when a > 10 deg . 

Note also that a(a-l0) l + and (<z-10)+ , for example, 

are in general not the same function. Because higher- 
order splines are created automatically by the 
multiplication of first-order splines (cf. Eq. (35)), only 
first-order splines in the explanatory variables need to 
be included as pseudo -explanatory variables. 

Extrapolating from this simple example, it is clear 
that multiple knots in the explanatory variables would 
provide a very effective nonlinear modeling capability 
in multiple dimensions with local nonlinear modeling 
capability, while retaining physical insight into the 
functional dependencies. The automated 

orthogonalization and sorting process described earlier 
identifies which nonlinear terms are necessary to 
characterize the functional dependencies and estimates 
(a) zero-degree, (b) first-degree, (c) second-degree the associated model parameters. Inputs required from 

the analyst relate only to the limits of what should be 
considered, such as which explanatory variables to consider, maximum model order to consider, and knot locations 
to consider. However, these can be specified very generously, because the orthogonal function modeling algorithm 
automatically sorts out which of the terms are important, based on the data, and discards the rest. The result is a 
global parsimonious model with excellent local nonlinear modeling capability and easy physical interpretation. 

IV. Test Aircraft and Flight Data 

A. T-2 Subscale Jet Transport Aircraft Description 

The T-2 aircraft is a 5.5 percent dynamically- scaled model 
of a generic commercial twin-engine jet transport aircraft. 

Figure 7 shows a photograph of the aircraft in flight. The 
aircraft has twin jet engines mounted under the wings and 
retractable tricycle landing gear. Aircraft geometry and 
nominal mass properties are given in Table 3. Further 
information on the T-2 subscale jet transport aircraft and 
associated flight test operations can be found in Refs. [24]-[26]. 

A similar airframe was tested extensively in the wind tunnel 8 , 
although the wind tunnel test article differed in some geometric 
details from the aircraft used for the flight tests. The wind 
tunnel data provided a reference for comparison with the 
aerodynamic models identified directly from flight data. 

7. Control Surfaces 

Control surfaces on the T-2 aircraft are left and right ailerons, left and right inboard and outboard elevators, 
upper and lower rudders, left and right inboard and outboard trailing-edge flaps, and left and right inboard and 
outboard spoilers, for a total of 16 independent control surfaces. For the flight data analyzed in this work, only the 
elevators, ailerons, and rudders were deflected. The individual elevator surfaces were moved together as a single 
elevator surface, and similarly for the rudders. Left and right ailerons were deflected asymmetrically, in the 



Figure 7. T-2 subscale jet transport aircraft 
Credit: NASA Langley Research Center 
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Figure 6. Polynomial splines in angle of attack: 
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conventional way. Definitions of control surface deflections are given below. Trailing edge down is positive 
deflection for wing and elevator surfaces, and trailing edge left is positive for rudder surfaces. 

< 41 > 

The aircraft can be flown by a safety pilot using direct visual contact and conventional radio control. A research 
pilot executed the flight test maneuvers from inside a mobile control room, using a synthetic vision display drawn 
from telemetry data and a local terrain database, along with video from a camera in the nose of the aircraft. Inputs 
from the research pilot and a ground-based flight control system were used to generate control surface commands 
which were transmitted by telemetry to the aircraft. 

The flight control system has the capability to inject automated control surface perturbations to excite the aircraft 
dynamic response for modeling purposes. These control surface perturbations can have arbitrary waveforms, and 
can be applied to multiple control surfaces individually or simultaneously. The perturbations are summed with pilot 
and feedback control commands in the flight control system, just before the actuator command rate and position 
limiting. Typically, the research pilot flies the aircraft to the desired flight condition(s), then initiates the automated 
control surface perturbations with a trigger switch on the control stick. 

2. Instrumentation and Data Acquisition 

The T-2 aircraft was equipped with a micro-INS, which provided 3-axis translational accelerometer 
measurements, angular rate measurements, estimated attitude angles, and GPS velocity and position. Air data 
probes attached to booms mounted on each wingtip (visible in Figure 7) measured angle of attack, sideslip angle, 
static pressure, and dynamic pressure. Measurements from static pressure sensors and ambient temperature sensors 
were used to compute air density and altitude. Engine speeds in rpm were measured and used as inputs to an engine 
model to compute thrust. The engine model was identified from ground test data, with adjustments for ram drag 
identified from flight data. Potentiometers on the rotation axes of the control surfaces measured control surface 
deflections. Mass properties were computed based on measured fuel flow, pre-flight weight and balance, and inertia 
measurements done on the ground for the aircraft without fuel. The pilot stick and rudder pedal commands and 
throttle position were also measured and recorded. Data from onboard sensors were telemetered to the ground in 
real time. Sampling rate for the flight data was 200 Hz, decimated to 50 Hz for data analysis and modeling. 

V. Simulation Results 

The techniques for flight test maneuver design and global aerodynamic modeling were applied to a nonlinear 
aircraft simulation of the T-2 subscale jet transport aircraft. The simulation was based on wind tunnel data 8 
collected using a wind tunnel model with geometry similar to the T-2 aircraft. Full nonlinear equations of motion 
were implemented and solved in MATLAB® using aircraft simulation software included in SIDPAC 14 , adapted for 
the T-2 aircraft. Documentation for the nonlinear aircraft simulation software can be found in Appendix D of 
Ref. [14]. The engine model was identified from ground test data, with ram drag corrections identified from T-2 
flight data. Simulated aircraft responses were corrupted with white Gaussian noise using root -mean- square 
amplitudes similar to what was observed in the T-2 flight data. 

Figure 8 shows control surface deflections and aircraft responses for a maneuver with orthogonal optimized 
multi-sine perturbation inputs applied to the elevator, aileron, and rudder control surfaces simultaneously, while the 
pilot implemented a slow increase in angle of attack by gradually pulling aft on the longitudinal stick with constant 
power setting. This is the flight test maneuver described in Section II for collecting data for global aerodynamic 
modeling. Figure 9 shows a cross-plot of angle of attack and sideslip angle, demonstrating how the angle of attack 
and sideslip angle are excited in an uncorrelated fashion, similar to Figure 4. Analogous cross -plots could be made 
for the other aircraft responses and control surface deflections. 
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Applying the global modeling techniques in simulation allows comparison of the global aerodynamic model 
identified from the simulated data with the known underlying aerodynamic database for the simulation. The 
aerodynamic database in the simulation was implemented in tabular form using wind tunnel data. 



time (s) time (s) time (s) 

Figure 8. Simulated flight data for a global aerodynamic modeling maneuver 


For the global aerodynamic modeling, initial modeling was done using coarse spacing of the knot locations in 
each explanatory variable to identify which explanatory variables required the increased local nonlinear modeling 
capability provided by spline terms. This analysis showed that the global modeling problem for the T-2 simulation 
data shown in Figure 8 required spline terms only for angle of attack. First-order splines in angle of attack with 

knots finely spaced at 6, 8, 10, 12, and 14 deg were 
added as pseudo-explanatory variables for the global 
aerodynamic modeling. The spline knot spacing 
could be made finer, with no penalty other than 
increased computation time to identify the model 
structure, as discussed earlier. For each 
nondimensional force and moment coefficient model, 
model terms with up to third-order complexity were 
orthogonalized and sorted to identify an adequate 
model structure and estimate the associated model 
parameters values and uncertainties. The third-order 
complexity applied to the conventional explanatory 
variables (such as angle of attack, nondimensional 
pitch rate, and elevator control deflection), and also to 
the first-order splines in angle of attack, which were 
included as pseudo-explanatory variables. 
Consequently, the first-order spline terms were also 
multiplied by themselves and each other, up to third 
order, so that this approach really allowed for up to 
third-order spline terms in the modeling. 

Figure 10 shows three-dimensional plots comparing the simulated flight data (x markers) with a mesh surface 
drawn by interrogating the aerodynamic database in the simulation and a smooth surface drawn using the global 
aerodynamic model identified from simulated flight data alone. The ranges of the explanatory variables used to 
generate the mesh surface and the smooth surface were chosen to include all values from the simulated flight data. 
The results in Figure 10 show that for all force and moment coefficients, the smooth surface and mesh surface are 
very similar, indicating that the global aerodynamic modeling approach captured the underlying functional 
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Figure 9. Simulated flight data cross plot for a 
global aerodynamic modeling maneuver 


dependencies accurately and efficiently, using data from a single flight test maneuver. This is true even in cases 
with complex local nonlinearity, as in the case of axial force coefficient C x and rolling moment coefficient Cj . 
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Figure 10. Global aerodynamic models identified from simulated flight data 


Close inspection shows only small mismatches between the smooth surface and mesh surface at locations 
relatively far from the simulated flight data, which would be areas of extrapolation for the identified global model. 
Flying another maneuver through that part of the explanatory variable subspace would provide the data necessary to 
correct those small mismatches. The three-dimensional plots shown in Figure 10 must necessarily omit some of the 
explanatory variables (because there are more than two explanatory variables, and the third dimension is used for the 
response variable), so the mismatch between simulated flight data and the smooth and mesh surfaces is mostly a 
dependence on an explanatory variable that cannot be shown in the three-dimensional plot. The results in Figure 10 
show that a few well-flown maneuvers of the type described here could very effectively cover a large portion of the 
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explanatory variable subspace and therefore could be used in place of wind tunnel testing that requires extensive 
grids of test points to generate the aerodynamic database. The approach described here could be extended to more 
than the three control surfaces used in this demonstration 16 , further increasing the productivity and efficiency of the 
flight test maneuver and the global aerodynamic modeling process. 

As an example of the final global modeling results, the global model identified for nondimensional aerodynamic 
vertical force coefficient C z from simulated flight data was 


C z ~ Cz. 


+ Cy CC + Cy — f Cy 

2 y $e 


5 e +Cz 


(, ■x-\0)\+C Zi (a -I2) l + 


(42) 


where model parameter estimates and uncertainties are given in Table 4. 

VI. Flight Test Results 

Flight test maneuvers similar to the one used for the simulation demonstration were flown on the T-2 subscale jet 
transport aircraft to collect data for global aerodynamic modeling. Figure 3 shows measured control surface 
deflections and aircraft responses for one of these maneuvers, with orthogonal optimized multi-sine perturbation 
inputs applied to the elevator, aileron, and rudder control surfaces simultaneously, while the pilot implemented a 
slow increase in angle of attack by gradually pulling aft on the longitudinal stick at constant idle power setting. The 
pilot also applied relatively large longitudinal and lateral stick inputs at roughly 20 sec, to recover the aircraft from a 
roll-off departure that occurred near 14 deg angle of attack. Automated perturbation inputs continued throughout the 
departure and recovery. 






Figure 11. Prediction using global aerodynamic models identified from T-2 flight data 


Detailed specifications of global aerodynamic models identified from flight data alone for all six nondimensional 
aerodynamic force and moment coefficients are given in Tables 5 and 6. These models were identified using the 
global aerodynamic modeling procedure applied to flight data from only the single maneuver shown in Figure 3. 
The identified global aerodynamic models were then installed in the nonlinear aircraft simulation for prediction 
testing. 

Figure 1 1 shows a prediction test using flight data from a piloted doublet sequence on longitudinal stick, rudder 
pedal, and lateral stick. The flight data used for this prediction test was of course not used in identifying the global 
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aerodynamic models. The solid lines show flight data and the dashed lines were computed by applying measured 
control surface deflections from the flight data to the nonlinear aircraft simulation, using the global aerodynamic 
models from Tables 5 and 6 for the aerodynamics. The results show reasonable accuracy in predicting the aircraft 
response for this maneuver, which would be classified as a local maneuver. This maneuver was chosen because it 
involved localized multi-axis aircraft response to elevator, rudder, and aileron deflections dissimilar from the inputs 
used to identify the global aerodynamic model. The predictions show that the global aerodynamic modeling 
approach applied to flight data from only a single flight test maneuver produced a reasonable prediction of the flight 
responses, indicating that the functional dependencies of the nondimensional force and moment coefficients on the 
explanatory variables were captured by the global aerodynamic modeling procedure. Presumably, using flight data 
from multiple runs of global aerodynamic modeling maneuvers similar to the one shown in Figure 3 would produce 
even better global aerodynamic modeling results. This is currently being investigated. In addition, it should be 
possible to identify thrust effects using global aerodynamic modeling maneuvers at different power settings, such as 
those listed in Table 2. 






Figure 12. Prediction using T-2 wind tunnel database 


Figure 12 shows prediction results for the same piloted doublet sequence maneuver, but instead using the wind 
tunnel database 8 for the aerodynamics in the nonlinear aircraft simulation. Figures 11 and 12 show similar 
prediction accuracy, suggesting that the global aerodynamic modeling procedure applied to data from a single flight 
test maneuver can produce simulation results with accuracy comparable to what can be obtained using aerodynamic 
data from an extensive wind tunnel test. In fact, the root-mean- square of the differences between the flight data and 
the predictions using the identified global model, shown in Figure 11, was 20 percent lower than for the predictions 
using the wind tunnel database, shown in Figure 12. Finally, the compact analytical form of the global aerodynamic 
model provides significant insight into the functional dependencies and is compact enough to be shown in Tables 5 
and 6. 
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VII. Concluding Remarks 

Simulation and flight test data for a subscale jet transport aircraft were used to demonstrate a novel technique for 
efficiently and effectively identifying global aerodynamic models from a single flight test maneuver. Novel flight 
test maneuver design, an orthogonal function modeling technique using splines, and an automated modeling process 
were combined to create an approach that very efficiently produced accurate global aerodynamic models with easy 
physical interpretation, based on flight data alone. The resulting global aerodynamic models were incorporated into 
a nonlinear simulation that exhibited reasonable prediction capability for data from flight maneuvers that were not 
used in the modeling process. 

The global modeling flight test maneuver was used to collect data over relatively large ranges of the explanatory 
variables for all six rigid-body degrees of freedom of the aircraft simultaneously. This resulted in highly efficient 
collection of flight data with de-correlated explanatory variables, for good model identifiability and modeling 
accuracy. The global aerodynamic modeling method incorporated spline functions in an automated orthogonal 
function modeling scheme, to allow accurate modeling of local complicated functional dependencies without 
adverse effects on other parts of the global model. 

Accurate global models identified from flight data alone can be used in rapid and relatively low-budget 
unmanned aircraft programs to save development time and money by making it possible to generate an accurate 
nonlinear aircraft simulation from a single flight, without the need for extensive wind tunnel testing or aerodynamic 
calculations. The capability also has important implications for aircraft safety, because the technique could be used 
to generate an onboard model of the global aerodynamics for an aircraft. The approach could be used to account for 
the particular geometry or flight environment of an individual aircraft, as well as provide a capability to monitor and 
account for changes in the aircraft due to failures, damage, and airframe icing, for example. A global aerodynamic 
model identified from flight data alone could also be compared to a simulation database generated from wind tunnel 
data and aerodynamic calculations for the purpose of improving the fidelity of these ground-based aerodynamic 
prediction methods. This could be done through identifying where there are significant differences between the 
database generated using ground-based methods and the global model identified from flight data alone, as well as 
where the models agreed well. Such comparisons could also be used to estimate full-scale Reynolds number effects, 
or artifacts of wind tunnel testing, such as sting interference and wall effects. 
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Table 1 Multiple input design for the T-2 subscale jet transport aircraft, 
level flight, V Q = 80 kts , a Q - 5 deg , T = 20 sec 


Input 

A (deg) 

\ (deg) 

k 

<fik ( rad ) 

RPF 

S e 

1.0 

0.3162 

5 

-2.2926 

1.13 

0.3162 

8 

0.6842 

0.3162 

11 

-0.3288 

0.3162 

14 

2.1677 

0.3162 

17 

-2.8795 

0.3162 

20 

-0.0447 

0.3162 

23 

-2.8485 

0.3162 

26 

-2.8634 

0.3162 

29 

3.0356 

0.3162 

32 

2.7574 

Sr 

2.0 

0.6325 

6 

0.9222 

1.04 

0.6325 

9 

0.7188 

0.6325 

12 

2.8103 

0.6325 

15 

-0.9035 

0.6325 

18 

0.1563 

0.6325 

21 

-1.8697 

0.6325 

24 

-0.8279 

0.6325 

27 

-2.0493 

0.6325 

30 

1.1970 

0.6325 

33 

0.5219 

S a 

1.0 

0.3162 

4 

1.8549 

1.17 

0.3162 

7 

-2.6561 

0.3162 

10 

-2.8832 

0.3162 

13 

0.1226 

0.3162 

16 

2.5070 

0.3162 

19 

2.6150 

0.3162 

22 

0.6119 

0.3162 

25 

-1.9709 

0.3162 

28 

-1.3854 

0.3162 

31 

-3.0152 
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Table 2. T-2 maneuvers applying multi-axis excitation for slow approach to stall and re covery 


Configuration 

Power Level 

Flight 

Card 

Cruise 

IDLE 

5 

10 

Cruise 

40 percent 

6 

11 

Cruise 

50 percent 

6 

12 

Cruise 

60 percent 

6 

13 

Takeoff 

IDLE 

6 

15 

Cruise 

IDLE 

9 

10 

Powered Approach 

IDLE 

9 

17 

Powered Approach 

IDLE 

11 

17 

Cruise 

IDLE 

34 

17 

Cruise 

IDLE 

34 

17 

Cruise 

IDLE 

36 

17 

Cruise 

IDLE 

36 

17 

Cruise 

IDLE 

37 

17 

Cruise 

IDLE 

37 

17 

Cruise 

IDLE 

38 

17 

Cruise 

IDLE 

38 

17 

Cruise 

IDLE 

38 

17 

Cruise 

IDLE 

42 

17 

Cruise 

IDLE 

42 

17 

Cruise 

IDLE 

42 

17 


Table 3. T-2 aircraft geometry and nominal mass properties 


c , ft 

0.915 

b , ft 

6.849 

S , ft 2 

5.902 

^,in 

57.30 

> in 

0.000 

z 0 ,in 

11.28 

, in 

56.63 

y cg > in 

0.000 

z cg , in 

11.43 

m , slugs 

1.585 

I x , slugs-ft 2 

1.179 

l y , SlUgS-ft 2 

4.520 

I z , slugs-ft 2 

5.527 

I xz , slugs-ft 2 

0.211 
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Table 4 . Parameter estimates and uncertainties for C z global model identified from simulated flight data 


C z Model 
Parameter 

Estimate 
± Std. Error 

C z 

-0.0738 

+0.0005 

C z 

-4.4809 

+0.0047 

C z 

^ q 

-47.241 

+0.2017 

Cz s e 

0.4500 


+0.0062 

c z 1 

2.1745 

“10 

+0.0124 

C Z 1 

0.8776 

a n 

+0.0167 
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Table 5. Parameter estimates and uncertainties for longitudinal global model identified from T-2 flight data 


C x Model 
Parameter 

Estimate 
+ Std. Error 

C z Model 
Parameter 

Estimate 
+ Std. Error 

C m Model 
Parameter 

Estimate 
+ Std. Error 

c x 0 

-0.1017 

+0.0039 

c z 

-0.1353 

+0.0037 

r 

m s e 

-0.8377 

+0.0822 

C *a 

1.6921 

+0.0754 

c z a 

-4.0013 

+0.0409 

r 

a 

6.5387 

+0.6531 

C x , 

a 2 

-4.4473 

+0.3301 

C Z 1 

*6 

1.0897 

+0.0488 

r 

v -"m i 
aag 

-1.7790 

+0.6789 

Cx q s e 

-306.56 

+15.289 

C z 1 

qa n 

2.0933e+03 

+134.16 

r 

m q 

-36.3742 

+4.9737 

C X ! 

qa lQ 

-806.31 

±53.991 

c 7 

A qS e 

214.2068 

+15.328 

r 

m G 

0.1592 

+0.0059 

c x s 

°e 

0.3982 

+0.0469 



r 

'm 

rn a 

-2.6947 

+0.1243 

c x 

^ aS e 

-1.7624 

+0.2961 



r 

m qS e 

388.9319 

+34.404 

C X ! 
a \o 

2.2999 

+0.1708 



r 

i 

qa n 

1.4879e+03 

+175.32 

c x , 

-7.9305 

+0.6953 



c 

m aq 

237.4375 

+50.250 





C m 

TYl i 

$e a \2 

2.7147 

+0.6377 





r 

i 

qag 

-198.14 

+73.932 





c 

m aS e 

-4.3485 

+0.7122 





r 

i 

S e ag 

5.5804 

+1.1554 
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Table 6. Parameter estimates and uncertainties for lateral global model identified from T-2 flight data 


C Y Model 
Parameter 

Estimate 
+ Std. Error 

Cj Model 
Parameter 

Estimate 
+ Std. Error 

C n Model 
Parameter 

Estimate 
+ Std. Error 

S 

-0.5483 

+0.0422 

C <* 

-0.0351 

+0.0009 

C np 

0.1673 

+0.0104 

C Y a 

0.0890 

+0.0051 

Cl rSa 

-0.6816 

+0.0690 

C n a 

-0.0017 

+0.0002 

Cy 

1 aS r 

0.5875 

+0.1537 

C lpr 

-2.3733 

+0.5028 


-0.1705 

+0.0020 

Cy s 

r5 a 

-3.7945 

+0.5046 

C lr 

0.2851 

+0.0282 

C “arS a 

5.6599 

+1.0124 

Cy 

± r 

2.4864 

+0.0528 

C l , 

aa io 

1.6212 

+0.0717 

C n 

n r 

-0.8890 

+0.0466 

Cy 

* pr 

6.3526 

+3.0915 

c v 

-0.1301 

+0.0046 

^ n apS a 

-1.3790 

+0.3426 

Cy 

P 5 a 

1.2747 

+0.1657 


5.6202e-03 

+0.0007 

c„ 

n a 

-0.0170 

+0.0013 

Cy 

± ap 

-10.1120 

+0.5171 

C 

-0.3172 

+0.0059 

C n 

n ap 

2.8815 

+0.1457 

Cy 

a5 a 

-2.3867 

±0.0928 

% 

0.0392 

+0.0011 

^ n ccS a 

0.3648 

+0.0221 

Cy 

‘fia 

-4.0360 

+0.2620 

c c 

a lO 

-0.3198 

+0.0164 

C n 

Up 

-0.4471 

+0.0275 

CY Sr 

0.2479 

+0.0256 

' ra lo 

3.2163 

+0.5092 

S, 

-0.0485 

+0.0031 

Cy ps a 

2.2167 

+0.1807 

C lpa 

-0.2710 

+0.0293 

^ n Pa 

0.3415 

+0.0665 

C Y S 

°a 

0.2640 

+0.0145 

l P a 10 

2.0470 

+0.1348 

Cn PaS a 

-3.1717 

+0.4148 

Cy 

L P 

1.2558 

+0.1045 

Z P“ 14 

-8.4945 

+0.7066 

C n 

ar 

1.9498 

+0.3031 

Cy 

5.3615e-03 

+0.0007 

Cl ps a 

0.1441 

+0.0176 

C,l pS r 

1.0237 

+0.1960 

Cy 

ps r 

-4.1527 

+0.7488 

Ci 

l o 

8.7560e-04 

+0.0001 





C i ar 

1.0919 

+0.2052 





C hs r 

1.3430 

+0.2878 
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